---
title: "Lab 2.2 Sample Solution"
author: "Student"
date: "11/13/2025"
format:
  pdf:
    colorlinks: true
---

```{r}
#| label: setup
#| include: false

# Load your dataset and place your libraries here
library(tidyverse)
library(broom)
library(kableExtra)
library(car)
library(sjPlot)
bicycle_data <- read.csv("seoul.bicycle.data.csv")
```

# Creating a multivariate model

For this assignment, we will be using data from the Seoul bicycle sharing system. In particular, we want to try and create a model of the `Rental.Count` of all bicycles rented per day.

> You should make a model that includes `Rental.Count` as the response variable and at least 2 quantitative variables as predictor variables. You can find the variable definitions [here](https://www.kaggle.com/datasets/joebeachcapital/seoul-bike-sharing).

## Forming expectations

1.  Before starting, write a sentence or two on your expectation of the direction and strength of the relationship between each of the predictor variables and the response variable

*I choose visibility and wind speed as the two predictor variables.*

*I expect a positive and strong relationship between visibility and rental count, since higher visibility usually indicates clearer weather and better riding conditions, encouraging more people to rent bicycles.*

*I also expect a negative and strong relationship between wind speed and rental count, as strong winds make cycling less comfortable and potentially unsafe, leading to fewer rentals.*

## Checking model fit

2.  Check the model fit via the `summary()` command information - does your model fit seem reasonable? Why or why not?

```{r}
#| label: run-regression
#| echo: false

model <- lm(Rental.Count ~ Visibility..10m. + Wind.speed..m.s., data=bicycle_data)
summary(model)
```

*The model generally matches my expectations: visibility shows a positive relationship with rental count, while wind speed shows a weak negative one. However, the R² value is very low (around 0.04), indicating that the model explains little of the variation in rentals. Overall, the fit is weak and not very reasonable.*

3.  Make a histogram of your residuals and a scatterplot of the residuals vs. fitted (or predicted). What does this indicate about your model fit? Does this seem like an appropriate model?

```{r}
#| label: fig-hist-residuals
#| fig-cap: "Histogram of residuals from Multiple Regression"
#| echo: false
#| message: false
#| warning: false

model.resids <- augment(model, bicycle_data)
ggplot(model.resids, aes(x = .resid)) +
  geom_histogram(fill = "lightblue", color = "black", alpha = 0.7) +
  labs(x = "Residual") +
  theme_minimal()
```

```{r}
#| label: fig-scatter-residuals
#| fig-cap: "Scatterplot of Residuals and Fitted Values"
#| echo: false
#| message: false
#| warning: false

ggplot(model.resids, aes(x = .fitted, y = .resid)) +
  geom_point(color = "blue", alpha = 0.7) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "Predicted values", y = "Residuals") +
  theme_minimal()
```

*The histogram of residuals looks bimodal rather than normally distributed, suggesting that the model doesn’t capture all underlying patterns in the data. However, the residuals vs. fitted plot shows no clear pattern, indicating that the variance is roughly constant. Overall, the model fit is only partially appropriate, it may miss some important predictors or nonlinear relationships.*

## Revising the model

4.  If the model does not appear to be a good fit, then revise it as best you can. If you have reached 30 min after lab started, go ahead and move on with whatever model you have so you will have time to complete the second half of the lab.

*Let's do a new models with Solar Radiation, Humidity percent, and Visibility.*

```{r}
#| label: run-regression2
#| echo: false

model2 <- lm(Rental.Count ~ Solar.Radiation..MJ.m2. + Humidity.percent. + Visibility..10m., data=bicycle_data)
summary(model2)
```

```{r}
#| label: fig-hist-residuals-new-model
#| fig-cap: "Histogram of residuals from New Model"
#| echo: false
#| message: false
#| warning: false

model2.resids <- augment(model2, bicycle_data)
ggplot(model2.resids, aes(x = .resid)) +
  geom_histogram(fill = "lightblue", color = "black", alpha = 0.7) +
  labs(x = "Residual") +
  theme_minimal()
```

```{r}
#| label: fig-scatter-residuals-new-model
#| fig-cap: "Scatterplot of Residuals vs. Fitted Values for the New Model"
#| echo: false
#| message: false
#| warning: false

ggplot(model2.resids, aes(x = .fitted, y = .resid)) +
  geom_point(color = "blue", alpha = 0.7) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "Predicted values", y = "Residuals") +
  theme_minimal()
```

*The histogram of residuals now looks approximately normal, suggesting that the model’s assumptions about normality are reasonable. The scatterplot shows a slight sign of heteroskedasticity, but the overall pattern appears fairly random, so the issue doesn’t seem severe.*

# Interpreting the coefficients

## Before we begin

1.  First, identify what are the units of the predictor and response variables?

*The response variable, `Rental.Count`, is measured in the number of bicycle rentals per day.*

*The predictor variables are:*

*Solar.Radiation..MJ.m2., measured in megajoules per square meter (MJ/m²),*

*Humidity.percent., measured in percent (%),*

*Visibility..10m., measured in units of 10 meters.*

2.  Second, make a table and identify what are reasonable values for the predictor variables and response variable. The table should include information about the 5 number summary for each.

```{r}
#| label: five-number-summary
#| fig-cap: "Five Number Summary for Key Variables"
#| echo: false

summary_table <- data.frame(
  Variable = c("Rental.Count", "Solar.Radiation..MJ.m2.", "Humidity.percent.", "Visibility..10m."),
  Min = c(min(bicycle_data$Rental.Count, na.rm = TRUE),
          min(bicycle_data$Solar.Radiation..MJ.m2., na.rm = TRUE),
          min(bicycle_data$Humidity.percent., na.rm = TRUE),
          min(bicycle_data$Visibility..10m., na.rm = TRUE)),
  Q1 = c(quantile(bicycle_data$Rental.Count, 0.25, na.rm = TRUE),
         quantile(bicycle_data$Solar.Radiation..MJ.m2., 0.25, na.rm = TRUE),
         quantile(bicycle_data$Humidity.percent., 0.25, na.rm = TRUE),
         quantile(bicycle_data$Visibility..10m., 0.25, na.rm = TRUE)),
  Median = c(median(bicycle_data$Rental.Count, na.rm = TRUE),
             median(bicycle_data$Solar.Radiation..MJ.m2., na.rm = TRUE),
             median(bicycle_data$Humidity.percent., na.rm = TRUE),
             median(bicycle_data$Visibility..10m., na.rm = TRUE)),
  Q3 = c(quantile(bicycle_data$Rental.Count, 0.75, na.rm = TRUE),
         quantile(bicycle_data$Solar.Radiation..MJ.m2., 0.75, na.rm = TRUE),
         quantile(bicycle_data$Humidity.percent., 0.75, na.rm = TRUE),
         quantile(bicycle_data$Visibility..10m., 0.75, na.rm = TRUE)),
  Max = c(max(bicycle_data$Rental.Count, na.rm = TRUE),
          max(bicycle_data$Solar.Radiation..MJ.m2., na.rm = TRUE),
          max(bicycle_data$Humidity.percent., na.rm = TRUE),
          max(bicycle_data$Visibility..10m., na.rm = TRUE))
)

summary_table %>%
  kbl() %>%
  kable_styling(full_width = FALSE)
```

3.  For what range(s) of our predictor and response variable(s) would our best fit line have meaning and when would it not? Can you think of any cases or examples of observation types that this model would not predict well?

*The model is only meaningful within the observed ranges of each variable. It wouldn’t predict well for extreme or unusual conditions beyond these ranges, such as abnormally high humidity or visibility values not seen in the dataset.*

## Interpreting the coefficients

3.  Interpret the intercept. Does it have meaning in this case?

*The intercept (-106.09) has no real meaning here, since it represents the predicted rental count when solar radiation, humidity, and visibility are all 0, conditions that do not occur in reality.*

4.  Interpret a one-unit change in each of the predictor variables on the response variable. Do you think this is a lot or a little? Why?

*Solar radiation: Holding other factors constant, a one-unit increase in solar radiation (1 MJ/m²) is associated with about 305 more rentals, which is quite a large effect.*

*Humidity: Holding others constant, a one-percent increase in humidity is linked to about 3.9 more rentals, which is a small effect.*

*Visibility: Holding others constant, a one-unit increase in visibility (10 m) adds roughly 0.04 rentals, which is negligible.*

5.  Make a table and solve for the predicted value ($\hat{y}$) of the response variable at Q1 and Q3 of each of the predictor variables, holding the other predictor variable(s) at their mean.

    | Predictions for $\hat{y}$ | $x_i$ at Q1 | $x_i$ at mean | $x_i$ at Q3 |
    |----------------------|-----------------|-----------------|-----------------|
    | Response var $x_1$ | `r -106.09335 + 305.20375 * 1.07 + 3.88782 * mean(bicycle_data$Humidity.percent.) + 0.04394 * mean(bicycle_data$Visibility..10m.)` | `r -106.09335 + 305.20375 * mean(bicycle_data$Solar.Radiation..MJ.m2.) + 3.88782 * mean(bicycle_data$Humidity.percent.) + 0.04394 * mean(bicycle_data$Visibility..10m.)` | `r -106.09335 + 305.20375 * 2.63 + 3.88782 * mean(bicycle_data$Humidity.percent.) + 0.04394 * mean(bicycle_data$Visibility..10m.)` |
    | Response var $x_2$ | `r -106.09335 + 305.20375 * mean(bicycle_data$Solar.Radiation..MJ.m2.) + 3.88782 * 34 + 0.04394 * mean(bicycle_data$Visibility..10m.)` | `r -106.09335 + 305.20375 * mean(bicycle_data$Solar.Radiation..MJ.m2.) + 3.88782 * mean(bicycle_data$Humidity.percent.) + 0.04394 * mean(bicycle_data$Visibility..10m.)` | `r -106.09335 + 305.20375 * mean(bicycle_data$Solar.Radiation..MJ.m2.) + 3.88782 * 57 + 0.04394 * mean(bicycle_data$Visibility..10m.)` |
    | Response var $x_3$ | `r -106.09335 + 305.20375 * mean(bicycle_data$Solar.Radiation..MJ.m2.) + 3.88782 * mean(bicycle_data$Humidity.percent.) + 0.04394 * 1083` | `r -106.09335 + 305.20375 * mean(bicycle_data$Solar.Radiation..MJ.m2.) + 3.88782 * mean(bicycle_data$Humidity.percent.) + 0.04394 * mean(bicycle_data$Visibility..10m.)` | `r -106.09335 + 305.20375 * mean(bicycle_data$Solar.Radiation..MJ.m2.) + 3.88782 * mean(bicycle_data$Humidity.percent.) + 0.04394 * 1948` |

    : Coefficient magnitude examination

6.  Does a change from the Q1 value to the Q3 of the predictor variable result in a large change (in your opinion) in the predicted value ($\hat{y}$) of the response variable?

*Based on the table, changing each predictor from its Q1 to Q3 value leads to a moderate to large increase in predicted rentals, especially for Solar Radiation​, which shows the largest change (about 477). This suggests that solar radiation has the strongest effect on rental counts compared to the other predictors.*

7.  Another way to view the same information is via a marginal effects plot. We can use the library `sjPlot` and the function `plot_model()` to assess how $\hat{y}$ changes as the predictor variable $x_i$ changes across its range, holding the other predictor variables at their mean value. You can view the documentation for the library [here](https://cran.r-project.org/web/packages/sjPlot/vignettes/plot_marginal_effects.html).

A sample command is the following:

```         
library(sjPlot)

mod <- lm(data=mtcars, mpg ~ wt + hp)

plot_model(mod, type="pred", terms="wt")
```

-   Create marginal plots of your predictor variables and interpret their meaning. Compare your output to the table you created in Steps 5 & 6.

```{r}
#| label: solar-rad-marginal-plot
#| fig-cap: "Effect of Solar Radiation on Rental Count"
#| echo: false
#| message: false
#| warning: false

plot_model(model2, type = "pred", terms = "Solar.Radiation..MJ.m2.")
```

```{r}
#| label: humidity-marginal-plot
#| fig-cap: "Effect of Humidity Percent on Rental Count"
#| echo: false
#| message: false
#| warning: false

plot_model(model2, type = "pred", terms = "Humidity.percent.")
```

```{r}
#| label: visibility-marginal-plot
#| fig-cap: "Effect of Visibility on Rental Count"
#| echo: false
#| message: false
#| warning: false

plot_model(model2, type = "pred", terms = "Visibility..10m.")
```

*The marginal plots show that rental counts increase strongly with solar radiation, moderately with humidity, and only slightly with visibility. These trends are consistent with the table results, solar radiation has the largest impact on rentals, while visibility shows the weakest effect.*

## Making a decision

8.  Decide whether each of your predictor variables has a substantively significant relationship with the response variable.

*Solar radiation has a substantively significant relationship with rental count, as higher sunlight levels are associated with increased rentals, and both the coefficient and marginal plot show a strong, meaningful effect.*

9.  If you have time, go back and try a more advanced model and repeat the above steps.
